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Abstract 

The analytic continuation of spectral functions is an important procedure to perform in several 
implementations of density functional theory (DFT) where many-body effects have been added 
through dynamical mean field theory (DMFT). Analytic continuation is numerically ill-posed, 
making this a non-trivial problem to solve. Here we investigate one of the most popular analytic 
continuation techniques, namely the Pade approximant. Aspects concerning its implementation 
are investigated with special regard towards making it stable and free of artificial defects. To this 
end electronic structure calculations are done using a single-channel scattering model, and the 
resulting DFT-level Green's functions are used to probe the properties of the Pade approximant. 
It is found that the analytic properties of the approximant can be controlled by appropriate 
modifications, making it a robust and reliable tool for electronic structure calculations. 
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I. INTRODUCTION 



Dynamical mean field theory combined with density functional theory formulated 
within the local density approximation (LDA+DMFT) is a rapidly developing area of the 
modern computational solid state physics of strongly correlated materials. A number of 
different implementations have been proposed and applied to calculate spectral functions 
and thermodynamic properties of Mott insulators, ferromagnetic 3d, 4/ metals, actinides 
or other systems Q]. Most of the LDA+DMFT calculations have been performed with 
partial self-consistency applied only to the many-body local self-energy, while the LDA 
potential (Hamiltonian) was kept fixed. In such schemes the effect of electronic correla- 
tions upon the charge is neglected. Recently, full self-consistent LDA+DMFT schemes 



became also available 



3j uses the full- 



The first implementation by Savrasov et al. 
potential linear muffin-tin orbitals (LMTO) method basis set, the many-body self-energy 
is added directly to the Hamiltonian, and the diagonalization of a non-Hermitian prob- 
lem is performed. Within the multiple scattering approach, the first implementation was 



realized within t 



: re exact muffin-tin orbitals (EMTO) basis set [4|, followed by the KKR 
implementation fj. While in the EMTO+DMFT approach the self-energy is added to the 
scattering path operator, from which the charge density is obtained by a complex contour 
integration, in the latter KKR+DMFT approach the self-energy is added directly to the 
radial Dirac solver from where the regular and irregular solutions are computed directly. 
In both these EMTO and KKR Green function schemes, the LDA solver is formulated 
within the complex plane. Therefore the combination with the many-body problem re- 
quires the "transfer" from the complex plane to the imaginary axis or the real axis where 
the many-body problem is solved. In both these implementations, the Pade analytical 
continuation is used to promote the Green's function from the complex contour into the 
many-body solver. 

In the present paper, we discuss several analytical continuation schemes, in order to 
clarify their accuracy in the LDA+DMFT calculations. We propose a new approach for 
an efficient and accurate scheme which is verified using the EMTO one-particle Green's 
function. Numerical tests are performed on solid hydrogen, which is a simple s-band 
system. 

The paper is organized as follows: in the next section we review the major implementa- 



tion steps of the EMTO and establish the one-particle Green's function that is subject to 
the analytical continuation. In Sec. IHIt we present details of the Pade analytical continu- 
ation and in Sees. II VlfVl we discuss its major implementation issues, as well as proposing 
procedures to test its accuracy. In Sec. [VI] we consider the example of hydrogen as a 
protoype s-band system to be investigated. 



II. OVERVIEW OF THE EMTO METHOD 



Within the multiple-scattering formalism, the one-electron Green's function G a (r, r', z) 
is defined for an arbitrary complex energy z as 



z + Vl- v\ ff (v)\ G a (r, r', z) = S(r - r'), (1) 

where a stands for the spin, and r and r' are the position coordinates. For most of the 
applications, e.g. standard KKR or LMTO methods, the LSDA effective potential from 
Eq. ([Q) is approximated by spherical muffin-tin (MT) wells centered at lattice sites R. 
Within a particular basis set, the one-electron Green's function is expressed in terms of 
the so-called scattering path operator, g^^z/W, as we ^ as the regular, Z RL (z,r R ), and 
irregular, J RL (z, t r ), solutions to the single site scattering problem for the cell potential 
at lattice site R, viz. 



g ^sda {tr + R rR/ + R , z) = j2 Z ° RL (z, r R )gt L 8 R,Uz)Z° RIL ,(z, r R> ) - 

L,L' 

- S RR/ ^r RL (z,T R )Z RL ( Z ,r R ), (2) 

L 

where L = (I, m) with / < l max (usually l max = 3) and = r — R denotes a point around 
site R. The real space representation for the scattering path operator for the muffin-tin 
potential is given by 

9Rl S £i>(z) = [5RR'hut a ~ RL {z) - Bbz^l^z)]- 1 , (3) 

where t RL (z) stands for the single scattering t- matrix and B R l^ R ili(z) are the elements of 
the so-called structure constant matrix. 



Conventional MT based KKR or LMTO methods have limited accuracy as a result of 
the shape approximation to the potential and charge density. The former method uses 
non-overlapping spherical muffin-tin potentials and constant potential in the interstitial, 
while the latter method approximates the system with overlapping atomic spheres ne- 
glecting completely the interstitial and the overlap between individual spheres. Recent 
progress in the field of muffin-tin orbitals theory p-llO] shows that the best possible rep- 
resentation of the full potential in terms of spherical wells may be obtained by using 
large overlapping muffin-tin wells 10Nl2| with properly treated overlaps. Within this 



so-called exact muffin-tin orbitals method 



6MlO| , the scattering path operator fl'L' ( z ) 



is calculated as the inverse of the kink matrix defined by 



k rl,r'L'( z ) = Srr'$lliD rl (z) - S RLtR , L i(z) 



(4) 



where D RL (z) denotes the EMTO logarithmic derivative function 

aa 

, and S RLjR <l>(z) is 

the slope matrix 6|. Since the energy derivative of the kink matrix, K RL R 'l'{z), gives the 
overlap matrix for the EMTO basis set 6|, the matrix elements of the properly normalized 



LSDA Green's function become 



9] 



.<-yct,LSDA 



'RL,R 



>L>{ Z ) 



9rl S rA''( z )^r''L'',r'L'( z ) - 5rr>Sll>Ir L ( z )i 



(5) 



R"L>< 



where I RL (z) accounts for the unphysical poles of K RLR , L ,(z). In the case of translation 
invariance, Eqs. (J4]) and ([5]), are transformed to the reciprocal space, so that the lattice 
index R runs over the atoms in the primitive cell only, and the slope matrix, the kink 
matrix, and the path operator depend on the Bloch wave vector k. 
The total number of states at the Fermi level Ep is obtained as 



N(E F ) = ^- E // G$£gb{k,z)dkdz, (6) 

Zm RL,R'L' J JBZ 

where the energy integral is carried out on a complex contour that cuts the real axis 
below the bottom of the valence band and at Ep. The k-integral is performed over the 
first Brillouin zone, which numerically is turned into a sum over a discrete k-mesh in the 
irreducible Brillouin zone (IBZ). 

The electronic structure calculation is started by an appropriate guess of the Fermi 
level and charge density, after which the semi-circular contour is created. A new Fermi 



level is calculated from Eq. fl6]), as well as a new charge density. From the charge density 
a new potential is constructed and this process is iterated until the Fermi level and the 
total energy of the system is converged below a set tolerance. 

In the EMTO+DMFT method, the many-body effects are added to the DFT-level 
Green's function through a local self-energy T, RL rl'( z ) v ^ a ^ ne Dyson equation 



Grl^R'L'O*-, z t 



GCT,LSDA 1-1 \ 
RL,R'L'\ K i z ) 



Srr'T, rlrl , 



(7) 



where G RL R , L ,(k., z) is now the LSDA+DMFT Green's function matrix. The local self- 
energy is obtained from a many-body solver defined on the fermionic Matsubara fre- 
quencies Uj = (2j + 1)ttT, where j = 0,±1,..., and T is the temperature, and it de- 
pends on the so-called effective medium or bath Green's function Q RL R , L i(ioS). This, in 
turn, is calculated from the k— integrated LSDA+DMFT Green's function, G RL R 'L'(ioo) = 
Ibz G RLiR/u {k,iu)dk, as 



RL,R'L' 



GRL,R'L'( iuJ ] 



^RR'^RL,RL'( iuJ ) 



Hence, the k— integrated LSDA+DMFT Green's function needs to be known on the Mat- 
subara frequencies. This is solved by analytically continuing the k— integrated Green's 
function from the semi-circular contour to the imaginary axis by the Pade method. In 
figure [H we illustrate the contours used in the calculations. In this work the focus will lie 
on the analytic continuation of the Green's function in the complex plane described in 
figure CD 



III. OVERVIEW OF PADE APPROXIMANTS 

An analytic function f(z) can be uniquely specified by an infinite number of its values 
at points contained in a compact set, on which the function is analytic. Construction of 
a sequence of rational functions f(z) having the same values as f(z) on a given set of 
points is known as the rational interpolation problem. 

The Pade approximant method is one possible way to solve the rational interpolation 
problem, and has proved very useful in providing quantitative information about the 
solution of many interesting problems in physics and chemistry. 

One of the most used techniques to construct the Pade approximant is the Thiele re- 
ciprocal difference method which was pioneered in the field of condensed matter physics 
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FIG. 1. Schematic picture of the complex energy contours used in the EMTO method. The 
charge density and number of states are calculated for points z on a semi-circle enclosing the 
valence band. For finite temperature calculations (e.g. LDA+DMFT) the equidistant Matsubara 
points iojj (in red) along the imaginary axis relative to the Fermi level Ep are used. The density 
of states is calculated for points z (blue) on a horizontal contour defined as K(z) + id (delta 
being a small real number). 



by Vidberg and Serene [l^] who used it as a means to analytically continue spectral 
functions from the imaginary axis to the real axis. The method of Thiele takes values 
f(zi),f(z2),...,f(zN) from A" (complex) points Zi,z 2 , Zn as input and returns an an- 
alytically continued function value f(z') at a point z'. This is performed with 0(N 2 ) 
complexity. If N is an odd number / will be a rational polynomial function where both 
the numerator and denominator will be of order (TV — l)/2. If N is even the denomintor 
polynomial will be of order N/2 while the numerator will be of order N/2 — 1. Hence, if 
the goal is to approximate a spectral function with M number of poles, then N = 2M 
input points should be used. In general M is not known a priori, and to resolve this issue 
Vidberg and Serene [ljj] recommended increasing iV until the approximant converges to 
a fixed value. 



Later on, the Pade approximant method was scrutinized by Beach et al. 15J who 
introduced the matrix formulation, which allows to compute the rational polynomial 
coefficients from which the approximant is constructed. As these coefficients are obtained 
by a matrix inversion the issue of accuracy is addresed using a multi-precision arithmetic 
in the calculation of the rational polynomial. In that study the approximant 

„Af-l 



/(*) 



ai + a 2 z + h a M z 



yM—l I Z M 



(9) 



b 1 + b 2 z + --- + b M z } 

was calculated explicitly by solving a linear system of (A^ even) equations to get the 
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complex polynomial coefficients {cii,bi}. In the following section we will present several 
issues concerning the use of the Pade method in electronic band structure calculations. 



IV. REVISED PADE APPROXIMATION 



A. k-integrated vs. k-resolved analytical continuation 



Pade 



G(k,z) 



+- C/(k, icj) 
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-*- G(lL0) 



FIG. 2. Schematic diagram of the analytic continuation before (straight path) or after (dashed 
path) Bloch summation. G(k, z) in the top left corner is given for a set of points. The goal is 



true function Q(iiS). Following the straight path, the Pade approximant needs to be constructed 
for each k 6 IBZ and the Bloch sum is performed at the desired point. Along the dashed path, 
only one single Pade approximant is needed. 

In the case of the k-resolved Green's function in Eq. [7J we have some a priori knowledge 
of the pole structure. This is because each orbital contributes with at least one pole 
(hybridization between bands will in general give rise to more than a single pole). On 
the other hand, the pole structure of the local Green's function in Eq. ([8]) is not as 
well defined. Since the Pade approximant method assumes a clear pole structure by 
construction, performing the analytic continuation on the k-resolved Green's function 
should be more suitable. The two different procedures are shown schematically in figure 



to obtain an analytic expression Q(iuj) (lower right corner) that can be used to approximate the 



El 
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B. Pole-zero pairs of the approximant 



An important issue is the nearly cancelling pole-zero pairs in the Pade approximant. 
These pairs arise when the order of the denominator in Eq. (Q is larger than the actual 
number of poles present in the function to be analytically continued. In the spirit of the 



matrix formulation of Beach et al. 



15] , in such cases the linear equations used to solve the 



rational interpolation problem becomes overcomplete. If this happens, the Pade algorithm 
tries to cancel the spurious poles by placing a zero of the numerator at the location of the 
pole. The pole-zero cancellation will not be exact due to the lack of arbitrary precision. 
These pole- zero pairs are defects in the approximant, since they are not part of the true 
function. Hence, these pairs should be found by performing a root-search on both the de- 
nominator and nominator polynomial and once found be removed from the approximant. 
Assuming that we have found a set of k pole- zero pairs {(pi, qi), (p2> Q2), •••> (Pfc> <Zfc)}> where 
the pj's are zeros and the q^s are poles of the approximant, the pairs can be filtered out 
by the following operation 

/»->/»n^- ( io ) 

f z-Pi 

In general it is difficult to make any clear statement on where exactly in the complex 
)lane these pole-zero pairs will appear. There are theorems about the overall distribution 



161 ] . but they are too general to be of any practical use for the present problem. 



1. Geometric search of true poles 

A pertinent question is how one is to decide which poles are physical and which are 
spurious. The most obvious way would be to make a geometric search in the complex plane 
and investigate whether a chosen pole has a zero in its neighbourhood. One of the simplest 
ways would be to let the neighbourhood be a circle of some predefined radius. The radius 
would have to be at least as large as the numerical accuracy used in the calculations, 
since this sets the minimum pole-zero distance. However, this straightforward method 
has some pitfalls. It is for example not obvious how to specify the neighbourhood in the 
most optimal way since the pole-zero pairs usually have different separations, as will be 
seen in Sec. IVII If the neighbourhood is defined too large then more than one zero can 
be present inside leading to ambiguity or the accidental removal of a physical zero. If on 
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FIG. 3. Schematic figure of a hypotetic geometric pole-zero search in the complex plane. The 
picture shows a true pole and a true zero of the approximant, located along the real axis. The 
rest of the poles and zeros are defects. To the left a too small neighbourhood (a circle with 
radius R') is choosen, not detecting the zero in the pole-zero pair. In the right part of the figure 
the neighbourhood (circle with radius R") is too large. This circle not only includes the zero 
belonging to the pole-zero pair, but also the zero coming from another pair as well as a true 
zero of the approximant. 



the other hand the neighbourhood is chosen too small, a spurious pole in a pole- zero pair 
with a large separtion could mistakenly be considered as a true pole of the approximant. 
Hence, in practice, this method is difficult to implement inside any self-consistent loop, 
where an algorithm would be necessary to automate the filtering of defects. A schematic 
picture illustrating the geometric search is shown in figure [3J 



2. Selection of true poles by multiplication with random numbers 

Another approach to select the true poles is to add a complex random number to the 
original input data jl^]. The reasoning behind this method is that the physical poles 
should be insensitive to this (numerical) perturbation while the spurious pole-zero pairs 
should redistribute themselves in the complex plane. Additive noise would effectively 
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lead to adding an extra zero to the nominator, destroying the asymptotic properties of 
the Green's function. In this study, we will chose to multiply the input Green's function 
with a real random number, in order to conserve the analytic properties of the function. 
Notice that multiplicative noise only leads to a rescaling of the pole weights, but keeps 
the analytic properties intact. Accordingly, the input Green's function is multiplied with 
real random numbers of the form 



.,N, 



where rj is a real scaling factor and the x^s are real random numbers between —0.5 and 
0.5. More detailed analysis of added noise in the approximant for simple functions can 



be found in Refs. 



18 



20]. 



3. The approximant as a sum of poles 

One can learn more about the Pade approximant by comparing it as a product of its 
poles and zeros, and as a sum of its poles and residues viz. 

M-l 

n (z-Pi) M 
K*) = c^ = E^t (12) 

i=l 

where Wi is the residue (or weight) of the z'th pole. C is a multiplicative complex constant 
that needs to be added back when reassembling the approximant as a product of its poles 
and zeros. This is because the root-finding procedure is invariant up to a scaling factor 
for the polynomial coefficients. The weights u>j can be found by some algebra. Multiply 
(JT2j) with (z - qj) 

M-l 

n {z-Pi) m 

C(z ~ Qj)^ = (Z- q 3 ) E — V (13) 

U(z-qi) qi 

i=i 

After some rearrangement we arrive at 

M-l 

n (z -p^ M 
n(z- qi ) & qi 
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If we now set z = qj we finally get that 

AI-l 

II (qj - Pi) 

"i = CJ W • ( 15 ) 

n (qj - Qi) 

Notice that all the poles have to be non-degenerate to avoid dividin g by zero. It is now 
possible to make a connection with two theorems (Theorem 43.1 in 2l|, and Ref. 221 ]) 
stating the relationship between a sum of the form in the right-hand side of Eq. (I12j) 
and its analytic properties. These theorems imply that if all the u>j's are positive and all 
the qi's are real, the approximant has positive-definite, integrable spectral weight and is 
analytic except along the real axis. This is what is needed for our spectral function to 
be physical. Using the above information, one can single out the physical poles from the 
defective ones, and write the approximant as 

/» = £'-^-, (16) 

i A I* 

where the prime indicates that the sum should be taken over physical poles only. Which 
poles that are true physical poles or defects should be decided on the ground of the above 
theorems. For example, if a pole has negative residue it can be considered unphysical and 
be removed. The Pade method will in general assign a non-zero imaginary part to the 
residues due to lack of arbitrary precision, even to the physical poles. Hence the exclusion 
of poles from the approximant due to non-zero imaginary parts has limited use. 



V. PROPOSED NEW SCHEME 



After consideration of the above issues, we propose and test the following scheme to 
perform the Pade analytic continuation of Green's functions in the EMTO method: 

• Calculate the k- resolved LSDA Green's function (the integrand in Eq. (jSJ)) and 
construct a Pade approximant for each k-point. Perform the Bloch sum only after 
the analytic continuation has been performed. 

• Search for physical poles by multiplication of a random number to the original input 
data and see how the poles and zeros redistribute themselves in the complex plane. 
In parallel, investigate the properties of all approximant pole weights. 

12 



The philosophy behind the proposed scheme is to perform analytic continuation on Green's 
functions with clearly defined poles. After the construction of the Pade approximant, the 
root-finding procedure is performed on the denominator and numerator polynomials using 
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231 ] . After the roots are found the pole- 



a Weierstrass-Durand-Kerner-Dochev method 
zero pairs are explicitly removed from the approximant. In the next section the above 
mentioned procedures will be investigated. 



VI. TESTS FOR THE ANALYTICAL CONTINUATION 

In the following section we will perform several numerical tests to validate the pro- 
posed method of the previous section. As a means to compare approximated and exactly 
calculated spectral functions we will consider the density of states as a function of energy 
DOS(E) = — ^Im[Tr(G(E + «0 + ))] of various systems. Using the Green's function cal- 
culated at the N first Matsubara frequencies, we will use the Pade method to analytically 
continue this function to points situated on a horizontal contour defined as z = 9ft(z) +i8, 
where 5 is a real constant (see figure [TJ). The Pade approximant will be constructed by 



the method used by Beach et al. [15| to obtain the rational polynomial form in Eq. fl9]). 
The test will be on solid hydrogen where the exact muffin-tin orbitals are expanded in a 
simple basis. We introduce an error Error, as a function of energy E 

Error(E) = \DOS{E) - DOS approx .(E)\ (17) 

as a measure of the deviation between the directly calculated density of states DOS(E) 
and the approximant DOS approx .(E) obatined from the Pade form. 



A. Solid hydrogen - a single s-band 

We will consider the simple system of hydrogen in the face-centered cubic structure. 
The EMTO basis will be a single s-state, avoiding any hybridization with higher ^-states. 
In the actual calculation, the Wigner-Seitz radius was set to 2.50 Bohr and the Bloch 
sum was performed using 10,569 k-points in the irreducible Brillouin zone. The EMTO 
equation was solved for 16 complex energy points along a semi-circular integration contour 
with a depth of 1.0 Ry relative to the Fermi level, enclosing the valence states. The slope 
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FIG. 4. Top: Density of states for solid hydrogen compared with approximants. The ap- 
proximants were constructed using the ^m-resolved local Green's function, using the N first 
Matsubara points corresponding to a temperature of 500 K. Inset: Error as defined in Eq. flTTJ). 
Bottom, left: Pole-zero distribution in the complex plane for the N = 16 approximant, giving 8 
poles and 7 zeros. 

Bottom, right: Pole-zero distribution in the complex plane for the N = 32 approximant, giving 
16 poles and 15 zeros. A pole at 2.0895 + i0.8561 Ry and a zero at 2.3022 + i0.2929 Ry not 
shown. 



matrix in Eq. was calculated using the two-center expansion [24[, and the local 
density approximation 25[ was used for the exchange-correlation functional. After self- 
consistency was reached the DOS was calculated as outlined above using a contour with 
8 = 0.02 Ry. 
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B. k-integrated Pade approximation 



First we investigate the Pade approximant for the k-integrated Green's function, i.e. 
the procedure illustrated by the dashed path in figure [2j The result obtained by ana- 
lytically continuing the local Green's function Gjn t n/Li(icUj), (L = L' = {£ = 0, m = 0}, 
R = R' single site), using the first N = 8, 16, 32 and 64 Matsubara points evaluated for a 
temperature T = 500 K are shown in figure H] (top panel). It can be noted that the Pade 
approximant can capture the smooth behavior of the DOS around the Fermi level, but it 
has difficulties in correctly reproducing the function at the band edges. Even though the 
approximant finds the correct location of the top of the band at roughly 0.1 Ry above 
the Fermi level, it either under- or overshoots the exact height while, at the same time, 
it introduces negative spectral weight into the DOS. 

C. Pade approximation for a single k-point 

Next we investigate the Pade approximant for a single k-point. For illustration pur- 
poses, here we select the T-point and test the Pade continuation for the Green's function 
calculated for k = (0, 0, 0). For the present test case (fee H), this spectral function has a 
physical pole at about -0.6 Ry below the Fermi level (corresponding to the bottom of the s 
band). As discussed in Sec. IIVP4 the defects in the Pade approximant should be sensitive 
to the addition of noise to the input data, while the physical poles should not. To test 
this effect, we multiply the input Green's function with random numbers as described in 
Sec. IIVBI The pole-zero distribution with and without multiplied random numbers are 
shown in the top of figure [5]for the N = 16 approximant using 77 = 10~ 6 . The positions of 
the poles and zeros from the original data are listed in table HI Except the (physical) pole 
(~ —0.6 Ry below the Fermi level) the 14 (unphysical) poles and zeros of the approxi- 
mant are scattered in the complex plane. 4 pairs are clustered around the imaginary axis, 
in the vicinity of the input points used to construct the Pade approximant. The other 
poles and zeros are situated further out in the complex plane. As the original input data 
is multiplied by random numbers, the distribution of the poles and zeros change. One 
extra pole- zero pair is now clustering around the imaginary axis in addition to the 4 pairs 
obtained for the original Green's function. All defects have changed their position while 
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the true pole has only been shifted slightly. That more of the pole-zero pairs wil 
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cluster 
26[, in 



around the vicinty of the input points was also observed by Sokolovski et al. 
that case for another kind of function. We should also mention that special care must 
be taken with the magnitude of the random numbers used in this procedure, as can be 
seen in the bottom panel of figure [5J In the bottom panel, where 77 = 10~ 3 , the true pole 
shows a larger shift then in the case of 77 = 10 -6 . Increasing 77 to even larger values than 
77 = 10 -3 gives a even larger shift of the true pole. 

Additional information about the approximant can be gained by investigating the 
residues of the poles, as explained in Sec. IIVB 31 The residues calculated according to 
(j!5p can be seen in table HJ together with the position of their respective poles. 4 of the 
residues Wi (i = 1,2,4 and 8) have negative real parts and their respective poles can be 
ruled out as spurious. Residues with label % = 3 and 6 have have real parts that are 
several orders of magnitude smaller than the residue i = 5, which has a real part close 
to 1. The residue of the i = 7 pole has a large imaginary part compared with the other 
poles, indicating that the residue violates the restrictions given by the theorems cited in 
Sec. IIVB 3[ Hence we should conclude that the pole q§ is our physical pole, and the rest 
are spurious. This conclusion is in perfect agreement with what was found above using 
random numbers on the input data. 



D. k-resolved Pade approximation 

Repeating the same procedure as the one for the k-integrated Green's function but 
this time for the k-resolved Grl^ijQs., iu)j), (L = L' = {£ = 0, m = 0}, R = R'), using 
N = 2, 4, 8 and 16 number of input points produces the result shown in figure O As seen, 
the N = 2 approximant fails to capture the DOS correctly at the bottom of the s-band. 
As more input points are added to the approximant a better fit is obtained. 

In figure [7] the DOS constructed using N = 2 approximants and varying the tempera- 
ture, T = 100, 250, 500 and 750 K is shown. This will change the position and distance 
between the input points used to construct the Pade. No qualitative difference is seen 
among the approximants, implying that it is not enough to go further from the real axis 
but that indeed a larger amount of input points are needed to approximate the parts of 
the DOS farther away from the Fermi level. 
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FIG. 5. Top: N = 16 approximant for the T-point, with and without noise. The true pole 
stemming from the original input data is situated at —0.5959 — i3.5719 • 10~ 6 Ry. Adding 
random noise with rj = 10 -6 gives the pole position at — 0.5963 + i3. 4499- 10 -5 Ry. The spurious 
poles and zeros show a larger scattering. The positions of the zeros and poles from the original 
data is tabulated in tabled! A zero of the true approximant at 4.8371 — i0.3177 Ry, and a zero 
from the noisy data at 3.2222 + z2.8388 Ry not shown. 

Bottom: Same as top, but here r] = 10~ 3 . The true pole has been shifted to —0.5803 — i2.0216 • 
10~ 2 Ry. 



Using the information given by the addition of random numbers and the residues, the 
approximant can now be reconstructed using only the true poles, see Eq. (fTBT) . We test 
this this reconstruction on the N = 16 approximant in figure El For each k-point the true 
pole was singled out. After this two options were tested. First, the sum in Eq. f fT6|) was 
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i Re{qi) 



Im(qi) 



Re(wi) 



Im(wi) 



f 1 1.1449469 • 10~ 2 0.1360169 -2.6324331 • 10" 13 -4.5109522 • 10~ 13 

f 2 -4.0764097 • 10" 2 0.3192908 -1.3977842 • 10~ 10 6.7675740 • 10~ n 

| 3 -0.5137978 1.111784 1.0184415 • 10~ 3 -1.3286583 • 10" 4 

f 4 4.0205983 • 10~ 3 0.2173124 -1.2889500 • 10" 13 -2.2621590 • 10" 13 

* 5 -0.5958811 -3.5719320 • 10" 6 0.9999138 4.7717123 • 10~ 5 

t 6 -2.8055273 • 10" 2 0.1243891 1.6093394 • 10" 12 -6.1991449 • 10~ 13 

t 7 1.547490 2.5824822 • 10~ 2 0.9696459 0.1238033 

f 8 2.441888 -0.1386907 -1.142889 -0.1035025 

Sum: 0.827689 2.021564 • 10" 2 

TABLE I. Real and imaginary parts of the N = 16 approximant poles with respective residues. 
The pole i = 5 (*) is a true pole of the approximant, while the rest (f) are spurious. The poles 
can be seen plotted in figure [5j 

kept as it is, with a fully complex weight. Second, only the real part of the weight Wi was 
used. In the second option we are enforcing the approximant to have the analytic structure 
specified by the theorems cited in Sec. IIVB 3[ which demands the residues to be purely 
real if our Green's function is to be physical. The relative error of both these procedures 
are plotted together with the original N = 16 approximant in figure El The reconstructed 
approximants shows a smaller error then the original approximant for energies above ~ 0.1 
Ry, and below ~ —0.6 Ry. It should also be noticed that the reconstructed approximant 
without any imaginary residue included in general gives a smaller error than the one 
where the imaginary part was kept. Looking closer around the Fermi level it can be seen 
that the reconstructed approximants are not always giving smaller errors then the original 
approximant. However, the point of the above procedure is not to improve the general 
fitting, but rather to ensure that no defects enter the approximant that would lead to 

18 




-1 -0.5 0.5 

E-E p [Ry] 



FIG. 6. Density of states for solid hydrogen compared with approximants. The approximants 
were constructed using the ^m-resolved k-dependent Green's function. Once the k-resolved 
Green's function was given for the horizontal density of states-contour the Bloch summation 
was performed. Inset: Error. 

non-analytic behavior. 



VII. SUMMARY AND CONCLUSION 

Pade approximants have long been in use for analytic continuation of spectral functions 
in condensed matter physics. This is since their rational polynomial form suits the ana- 
lytical properties of the spectral functions of interest. The instabilities often encountered 
when using these methods should be attributed to the ill-posedness of analytic continu- 
ation in general, and not to the Pade method itself. Indeed, as has been shown earlier 
[jjj] the method can reach remarkable accuracy when the conditions concerning input 
precision and arithmetics are good enough. Typically this is not the case for electronic 
band structure calculations. However, this study shows that the analytic continuation 
can give acceptable results if special consideration is taken. The biggest advantage given 
is that calculations are performed away from the real axis, making the exact spectral pole 
structure unnecessary to know. The important issue to resolve is to be able to ensure 
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FIG. 7. Density of states for solid hydrogen compared with N = 2 approximants using temper- 
atures T = 100, 250, 500 and 750 K. The approximants were constructed using the ^m-resolved 
k-dependent Green's function. Inset: Error. 

that no spurious poles or zeros enter the half-plane in which the contour integrations and 
many-body equations are solved. As shown above this can be done once the locations of 
the approximant poles and zeros are known, as one is then able to use physical reasoning 
and the properties of the approximant to remove any unwanted features. 
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FIG. 8. Density of states and the N = 16 approximant (green) compared with approximants of 
the same order, but reconstructed using Eq. (|16p . For one approximant (blue) the imaginary 
part of the residue was kept, for the other (violet) only the real part was used. Inset: Error. 

tional rescources. 
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